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We propose a selfconsistent microscopic model of verti- 
cal sequential tunneling through a multi-quantum well. The 
model includes a detailed description of the contacts, uses 
the Transfer Hamiltonian for expressions of the current and 
it treats the Coulomb interaction within a mean field approx- 
imation. We analyze the current density through a double 
well and a superlattice and study the formation of electric 
field domains and multistability coming from the Coulomb 
interaction. Phase diagrams of parameter regions (bias, dop- 
ing in the heterostructure and in the contacts, etc) where the 
different solutions exist are given. 
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Coulomb interaction in heterostructures with large 
area wells is a small effect compared with the energy dif- 
ference between non-interacting eigenstates of the struc- 
ture. Therefore a mean field model gives, for many pur- 
poses a good description of the system. Among fea- 
tures of the transport properties having their origin in 
Coulomb interaction, intrinsic bistability has great im- 
portance. This physical phenomenon arises from the non- 
linear effect of the electric charge on the induced elec- 
trostatic potential, and it has been predicted and ob- 
served in double barrier structures (DB) Further- 
more in the presence of a laser polarized in the sample 
growth direction, new bistability regions have been the- 
oretically predicted ||. In this paper we deal with the 
statics and dynamics of vertical transport through bi- 
ased heterostructures whose main mechanism is sequen- 
tial tunneling. This is a topic which has attracted a great 
deal of attention in recent times. In weakly coupled su- 
perlattices, multistability due to domain formation has 
been much studied both theoretically and experimentally, 
|^-^|. When the charge in the superlattice is small due 
to lower doping in the wells, self-sustained current oscil- 
lations and chaos due to domain dynamics are possible 
JlO| pjj . So far, the most successful modeling of these 
phenomena use discrete rate equations for the electron 
density and electric field in each well, plus constitutive 
laws for the current, bias, boundary and initial condi- 
tions, fl||Jl3[- The laws may be phenomenological || or 
obtained from microscopic considerations, [f7| ]l^|l5fl . In 
all cases cited, the boundary conditions were selected in 



a more or less ad hoc manner by using the available infor- 
mation from experiments. This is particularly annoying 
because the boundary conditions select the relevant dy- 
namics of electric field domains in the oscillatory regime 
|]l6f . In this paper we present a microscopic model which 
includes in a natural way boundary conditions due to 
the emitter and collector regions of a multiwell structure 
(MW). We then solve it for the cases of a double quan- 
tum well (DQW) and a superlattice (SL). The presence 
of intrinsic bistability is demonstrated through phase dia- 
grams and I-V characteristics obtained by numerical sim- 
ulation and by means of numerical continuation of sta- 
tionary solution branches. The main ingredients of our 
model are as follows: we assume that the characteristic 
time of intersubband relaxation due to scattering (about 
0.1 ps for optical phonon scattering Jl7t) is much smaller 
than the tunneling time (less than 0.5 ns), which is in 
turn much smaller than the dielectric relaxation times re- 
sponsible for reaching a steady state (about 10 ns for the 
9nm/4nm GaAs/AlAs superlattices of Ref. @). This 
separation of time scales, as well as the configuration 
of a typical sample allows us to consider that only the 
ground state of each well is populated and that the tun- 
neling processes are stationary. Our assumptions then 
justify using rate equations for the electron densities at 
each well with relations for the currents calculated by 
means of the Transfer Hamiltonian (TH) |ll| . The rate 
equations for the electron densities imply that the inter- 
well currents and the currents from the emitter and to 
the collector are all equal to the total current in the sta- 
tionary case (a form of Ampere's law may be derived in 
the time-dependent case). Furthermore, since no current 
is created or destroyed in the MW, the total charge in it 
(emitter and collector included) is zero. Finally, the elec- 
trostatics is simplified by assuming that the charges are 
concentrated on 2D planes located at the wells, emitter 
and collector regions, as indicated in Fig. [l] and further 
explained below. 

The hamiltonian for independent electrons in a MW 
under dc bias is: 
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Here ^.(c/^) are the creation (annihilation) operators in 

the leads and dWd^) are the creation (annihilation) op- 
erators in the wells, and T^. are the tunneling matrix 
elements. The latter depend on the local electric field and 
must be calculated self-consistently for each bias. Apply- 
ing the TH under the assumptions listed before, we ob- 
tain the following expressions for the tunneling currents 
where J e ,i and Jn,c are the currents in the contacts and 
Ji,i+\ the interwells currents: 
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where i — 1, . . . , N — 1, iV is the number of wells, n is the 
number of subbands in each well i with energies € % C a (ref- 
ered with respect to the origin of potential drops), and 
Tj are the transmision coefficients through the ith bar- 
rier. The spectral functions of the wells are Lorentzians 
whose widths correspond to the LO phonon lifetimes (~ 
1-10 meV): A 4 Cj (e) = 7 /[(e - e^) 2 + 7 2 ] for the ith well. 
Of course this model can be improved by calculating mi- 
croscopically the self-energies, which could include other 
scattering mechanisms (e.g. interface roughness, impurity 
effects |n|) or even exchange-correlation effects (which 
affect the electron lifetime in a self-consistent way [pS[). 
We have assumed that the electrons in each well are in 
local equilibrium with Fermi energies e Ui which define the 
electronic densities n,. For a given set {e Wi } the densities 
evolve according to the following rate equations: 
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In these equations J M+ i = Ji,i+i(e Ui , e m+1 , $), J e> i = 
Je,i(tu)i j 3>) j and Jn, c = Jn,c{^lj n ,^), where $ denotes 



the set of voltage drops through the structure which are 
calculated as follows. The Poisson equation yields the 
potential drops in the barriers, Vi, and the wells, V W i 
(see Fig. |): 
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where e is the GaAs static permittivity, ni(e Ui ) is the 2D 
(areal) charge density at the ith well (to be determined), 
w and d are the well and barrier thickness respectively, 
and N% is the 2D intentional doping at the wells. The 
emitter and collector layers can be described by the fol- 
lowing equations || : 
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To write the emitter equations (^|) , we assume that there 
are no charges in the emitter barrier. Then the electric 
field across Si (see Fig. |l|) is equal to that in the emitter 
barrier. Furthermore, the areal charge density required 
to create this electric field is provided by the emitter. 
N(Ep) is the density of states at the emitter E F - To 
write the collector equations (^), we assume that the re- 
gion of length 82 in the collector is completely depleted 
of electrons || and local charge neutrality in the region 
of length (53 between the end of the depletion layer 82 
and the collector. In order to close the set of equations 
we impose global charge conservation and that all volt- 
age drops across the different regions must add up to the 
applied bias: 
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Note that the right hand side of Eq.(||) is the positive 
2D charge density depleted in the collector region. In- 
stead of the rate equations (^) , we can derive a form of 
Ampere's law which explicitly contains the total current 
density J(t). We differentiate (||) with respect to time 
and eliminate i%i by using (0). The result is 



£ dl£ +Ji - 1 ' i - J{t) > 



l,...,N + \, (10) 



where J(t) is the sum of displacement and tunneling 
currents. The time-dependent model consists of the 
3N + 8 equations - (^fj) (the currents are given by 
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Eqs. (0)), which contain the 3iV + 8 unknowns e u i, V w 



,N), Vj (j = 1,...,N + 1), Ax, A 2 , 4 



(k = 1, 2, 3), er, and J. Thus we have a system of equa- 
tions which, together with appropriate initial conditions, 
determine completely and self-consistently our problem. 
The boundary conditions arise in a natural way. Notice 
that the charge and electric field at the boundaries can- 
not be set prior to the calculation of the whole structure, 
which all previous models did @§|n§. 
In this paper we are interested in analyzing the statics of 
the model and the stability of the stationary solutions. 
One way to do this is to numerically solve the algebraic- 
differential system (Q) - ([To]) (plus appropriate initial con- 
ditions) for each bias until a stationary profile is reached. 
This is rather costly, so that we follow this procedure for 
a given value of the bias and then use a numerical contin- 
uation method to obtain all stationary solution branches 
in the I-V characteristic diagram. This yields both un- 
stable and stable solution branches. Direct integration 
of the stationary equations [dropping the displacement 
current in (|l0|)1 presents important problems of numeri- 
cal convergence to the appropriate solutions in regions of 
multistability (see below). 

We analyze a DQW (sample a) consisting of 90A GaAs 
wells and AOA Ga.5Al.5As barriers. The doping at both 
emitter and collector is No = 2 x 10 18 cm~ 3 , and in the 
wells it is N% = 1.5 x 10 n cm- 2 . The half-width of the 
well states is 7 = 4meV in Eqs. (Q) and T = 0. We do 
not consider the effect of other symmetry points in the 
conduction band than T. Fig. |^ shows the DQW I-V 
characteristic for two different values of N%. In Fig. |2|a 
the low bias peak corresponds to C1C1 tunneling (Ci 
are the conduction subbands ordered starting from that 
with lowest energy) between adjacent wells. At higher 
bias multistability of stationary solution branches sets in 
(three stable solutions coexist at about 0.44 V). To un- 
derstand the difference between these solutions, we have 
depicted in Fig. || the potential profile of three different 
solutions (two stable, one unstable) corresponding to the 
same voltage (0.41 V). In Fig. ||a, the electrons flow from 
the emitter to CI in the first well. Then there is C1C2 
tunneling to the second well. A stable solution with lower 
current density is shown in Fig. |^c. There the CI at 
the first well is below the bottom of the emitter layer, 
then the electrons flow to C2 at the first well instead and 
J is smaller. A similar situation occurs at higher bias, 
Figs. ||d to f. The subbands of the two wells are clearly 
off-resonance for the solution with lowest current, Fig. |^f. 
Notice that the current flowing from the emitter to the 
structure and the potential profile are quite different for 
the different solutions of our model. This shows that 
boundary conditions assumed in previous publications 
might constitute gross oversimplifications of the physical 
situation. In Fig. ^ we present the regions of multistabil- 
ity depending on the bias and the dimensionless doping 
inside the wells, v = ew N^/[e{ec2 — eci)], or at the 



emitter s = ew 2 Njj /[e(ec2 — e ci)]- We see that there is 
a lower and upper limit for both v and s to have bistabil- 
ity. Then it is possible to control the presence and extent 
of bistability in a sample by changing the doping in the 
wells or at the contacts and the well widths. 

In Fig. I we plot the I-V 
curve of a 90AGaAs/40AGa. 5 Al. 5 As SL with 11 barri- 
ers and 10 wells. Its doping is as in Fig. ||a. The stable 
branches are shown as continuous lines in Fig. ||. The 
inset shows three electric field profiles corresponding to 
three different voltages. They show the presence of do- 
mains in the SL with a domain wall which moves one well 
as we change the bias from one branch to the next one. 
Domain coexistence is also shown in the SL electrostatic 
potential profile; see Figjj] for a fixed bias V2 = 0.81V. 
The first branch in Fig. ^corresponds to C1C1 tunnel- 
ing. As V increases, C1C2 tunneling becomes possible 
in part of the structure and we have domain formation. 
This situation confirms the findings with other discrete 
models with ad hoc boundary conditions |],[l3] 15 1. An 
interesting feature in Fig. || is that satellite peaks have 
a smaller current than the C1C1 peak. This agrees with 
the results of Ref. |15|. Another interesting feature due 
to the voltage drop at the contacts is that the number 
of branches in the I-V curve is less than the number of 
wells. This behaviour can be understood looking at the 
branch at 1.21 V where the low field domain occupies 
the two wells closer to the emitter. C1C2 tunneling oc- 
curs between all the wells in the branch with V3 = 1.48V 
corresponding to a very intense peak of the current. 

In summary, we have proposed and solved a micro- 
scopic selfconsistent model for the sequential current 
through a multiwell structure which includes the cur- 
rent through the contacts and appropriate boundary con- 
ditions. We have obtained the static I-V curve and 
phase diagrams of a DQW and a SL, which display mul- 
tistability associated to domain formation. Exchange- 
correlation (not included in our model) has been demon- 
strated to reduce the bistability in a DB |l| . Including 
exchange-correlation effects is the aim of a future work. 
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FIG. 1. Electrostatic potential profile in a SL (sample b) 
for V 2 = 0.81V. 

FIG. 2. DQW I-V characteristics (sample a). The 
continuous (dotted) lines correspond to stable (unstable) 
solution branches: (a) — 1.5 x 10 n cm~ 2 ; and (b) 

N% = 4.31 x 10 n cm~ 2 . The inset magnifies the CTCT reso- 
nant peak, showing the region of bistability. 

FIG. 3. (a) - (c) The three stationary potential profiles for 
the DQW structure of Fig. 2a at 0.41 V, ordered from highest 
to lowest current density, (d) - (f) Same for a bias V = 0.56V. 
In all cases the emitter Fermi energy and the subband energies 
are depicted. 

FIG. 4. Phase diagrams showing the regions of multista- 
bility for the DQW (sample a): (a) dimensionless well doping 
v versus voltage at si = 1.97 (N D = 2x 10 18 cm~ 3 ),^i = 0.46 
corresponds to the well doping of Fig. 2b; (b) dimensionless 
contact doping s versus voltage, at v\. 

FIG. 5. I-V characteristic curve of a SL (sample b). The 
inset shows the electric field distribution through the SL for 
three voltages: Vi = .69 V; V 2 = .81V; V 3 = 1.48V. 



4 




0.0 0.5 v 2 1.0 1.5 

Voltage(V) 



